BZOJ3231【SDOI2008】递归数列 <矩阵快速幂>

Problem

【SDOI2008】递归数列

Description

一个由自然数组成的数列按下式定义：

对于：
对于：

其中和是给定的自然数。写一个程序，给定自然数, 计算, 并输出它除以给定自然数的余数的值。

Input

输入由四行组成。
第一行是一个自然数。
第二行包含个自然数。
第三行包含个自然数。
第四行包含三个自然数。

Output

输出一行一个正整数，表示的值。

Sample Input

2
1 1
1 1
2 10 1000003

Sample Output

HINT

对于的测试数据：，

标签：矩阵快速幂

Solution

简单矩阵快速幂优化递推。

将答案拆成两个前缀和，即。我们需要快速计算的值，显然可以构造递推矩阵。

特判暴力，其余用转移矩阵的次方乘由和组成的矩阵即可得到。

Code

#include <bits/stdc++.h>
#define MAX_N 20
using namespace std;
typedef long long lnt;
template <class T> inline void read(T &x) {
    x = 0; int c = getchar(), f = 1;
    for (; !isdigit(c); c = getchar()) if (c == 45) f = -1;
    for (; isdigit(c); c = getchar()) (x *= 10) += f*(c-'0');
}
int n; lnt l, r, b[MAX_N+5], c[MAX_N+5], s1, s2, P;
struct Matrix {
    int n; lnt ele[MAX_N][MAX_N];
    Matrix (int _n) {n = _n, memset(ele, 0, sizeof ele);}
    inline Matrix operator * (const Matrix &x) const {
        Matrix ret(n);
        for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) for (int k = 0; k < n; k++)
            ret.ele[i][j] = (ret.ele[i][j]+ele[i][k]*x.ele[k][j]%P)%P;
        return ret;
    }
} ;
Matrix Power(Matrix mat, lnt k) {
    Matrix ret(mat.n);
    for (int i = 0; i < ret.n; i++)
        ret.ele[i][i] = 1;
    for (; k; k >>= 1, mat = mat*mat)
        if (k&1) ret = ret*mat;
    return ret;
}
int main() {
    read(n); Matrix p(n+1), q(n+1);
    for (int i = 1; i <= n; i++) read(b[i]);
    for (int i = 1; i <= n; i++) read(c[i]);
    read(l), read(r), read(P), l--, p.ele[0][0] = 1;
    for (int i = 1; i < n; i++) p.ele[i+1][i] = 1;
    for (int i = 1; i <= n; i++) p.ele[0][i] = c[i]%P;
    for (int i = 1; i <= n; i++) p.ele[1][i] = c[i]%P;
    for (int i = 1; i <= n; i++) (q.ele[0][0] += b[i]) %= P;
    for (int i = 1; i <= n; i++) (q.ele[n-i+1][0] = b[i]) %= P;
    if (l <= n) for (int i = 1; i <= l; i++) (s1 += b[i]) %= P;
    else s1 = (Power(p, l-n)*q).ele[0][0];
    if (r <= n) for (int i = 1; i <= r; i++) (s2 += b[i]) %= P;
    else s2 = (Power(p, r-n)*q).ele[0][0];
    return printf("%lld\n", (s2-s1+P)%P), 0;
}